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ARTICLEINFO ABSTRACT 


Keywords: The new emerged infectious disease that is known the coronavirus disease (COVID-19), which is a high conta- 

Fractional Caputo derivative gious viral infection that started in December 2019 in China city Wuhan and spread very fast to the rest of the 

oe 2i. . world. This infection caused million of infected cases globally and still pose an alarming situation for human 

Non-pharmaceutical intervention ; I : p ae ; , R ] š ; 

Stability lives. Pakistan in Asian countries is considered the third country with higher number of cases of coronavirus with 
abit 


more than 200,000. Recently, many mathematical models have been considered to better understand the 
coronavirus infection. Most of these models are based on classical integer-order derivative which can not capture 
the fading memory and crossover behavior found in many biological phenomena. Therefore, we study the 
coronavirus disease in this paper by exploring the dynamics of COVID-19 infection using the non-integer Caputo 
derivative. In the absence of vaccine or therapy, the role of non-pharmaceutical interventions (NPIs) is examined 
on the dynamics of theCOVID-19 outbreak in Pakistan. First, we construct the model in integer sense and then 
apply the fractional operator to have a generalized model. The generalized model is then used to present the 
detailed theoretical results. We investigate the stability of the model for the case of fractional model using a 
nonlinear fractional Lyapunov function of Goh-Voltera type. Furthermore, we estimate the values of parameters 
with the help of least square curve fitting tool for the COVID-19 data recorded in Pakistan since March 1 till June 
30, 2020 and show that our considered model give an accurate prediction to the real COVID-19 statistical cases. 
Finally, numerical simulations are presented using estimated parameters for various values of the fractional order 
of the Caputo derivative. From the simulation results it is found that the fractional order provides more insights 
about the disease dynamics. 


Parameter estimation 
Numerical simulations 


1. Introduction 


A highly contagious viral disease caused by a coronavirus (SARS- 
CoV-2), first identified in China at the end of 2019 and affecting 213 
countries and territories around the world [22]. The novel Wuhan virus 
figures under different names such as the novel coronavirus (2019- 
nCoV) on January 7 and COVID-19 by the World Health Organization 
(WHO) on February 11 [38]. The transmission of the virus was consid- 
ered initially animal-to-human and fist human-to-human transmission 
confirmed on January 20 in Guangdong, a coastal province of southeast 
China [38,12]. Over 10.27 million confirmed cases and more than 0.5 
million cumulative deaths all over the world are recorded till June 30, 
2020 [1]. The mortality rate is about 796 and the recovery rate is 9396 
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globally [1,22]. The pharmaceutical companies, regulatory authorities 
are doing every thing possible for the availability of a safe and effective 
anti-viral vaccine against this novel infection. The pandemic continu- 
ously inflects severe public health, hit the laborer’s community and 
economy throughout the world, including in Pakistan. The suggested 
incubation period for this viral infection is in the range 2 to 10 days by 
WHO [25,12]. The symptoms of a COVID-19 infected person include 
fever, myalgia or fatigue, dry cough, shortness of breath or dyspnoea, 
pneumonia, chills, sore throat, multiple organ failure, anorexia or lungs 
failure, acute respiratory distress syndrome (ARDS), lymphopenia (a 
reduction of lymphocytes in the circulating blood), loss of smell or taste 
and a sign of cytokine storm [25,12]. The majority of deaths have 
occurred in other chronic illnesses like diabetes, hypertension and 
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cardiovascular disease. In the absence of therapies, approved vaccine 
and antiviral the whole world is focusing on non-pharmaceutical in- 
terventions to eradicate this COVID-19. The most common NPIs prac- 
tices these days are social physical distance, exposed to be a quarantine, 
isolation of infected, wearing mask, more testing facility, closure of the 
school, non-essential business, hand washing or sanitizing, lockdown 
(stay at home and work from home) and protective kit for medical 
personal and avoiding unnecessary gathering. 

Pakistan is one of the epicenters of COVID-19 in Asia with an esti- 
mated population of over 220 million [1]. Currently, Pakistan is the 3rd 
country in Asia and 12th in the world with high confirmed infected cases 
[17,1]. As on 30 June 2020, 209,337 total confirmed cases are recorded 
and about 4304 lost their lives. Although in Pakistan the mortality rate is 
very low as compared to the recovery rate, still this disease has a 
negative impact on the economy. The first case reported in Karachi on 
26 Feb 2020, the virus spread quickly within three weeks in the entire 
country. Also, the accounted cases are less than the total cases due to 
limiting testing. Currently, the situation is worse, due to open lockdown 
on 9 May, the government is unable to maintain strict lockdown because 
of severe economic hardship. But the government imposed smart lock- 
down and easing restrictions, by allowing offices, business, markets with 
limited hours and staff capacity, five days in a week to reopen the 
economy of the country. 

Mathematical models coupled with statistical data and its analysis 
are better to obtained reasonable information about the coronavirus 
disease and its complications and future spread and control. Recently, it 
is a huge challenge that without vaccines, how to curtail the pandemic 
with non-pharmaceutical intervention (NPIs). In this regard, different 
mathematical modeling approaches have been adopted to analyze the 
transmission pattern of COVID-19. Ferguson et al. [16] studied a COVID- 
19-induced mortality model with the impact of NPIs. A mathematical 
model for accessing the community-wide impact of mask in controlling 
the spread of coronavirus is developed by Eikenberry et al. [14]. The 
dynamics and mitigation strategies for COVID-19 in Canada have been 
studied in [39]. A compartmental system describing the coronavirus 
dynamics in three highly infected countries has been suggested in [15]. 
Recently, Ullah and Khan studied the coronavirus model for Pakistani 
data using optimal control theory [32]. Fractional calculus is the 
generalization of classical calculus. Mathematical models with non- 
integer order operators provide a better understanding of a phenom- 
ena. Further, models with fractional order derivatives are capable to 
capture the fading memory and crossover behavior and provides a 
greater degree of accuracy. Mathematical models with fractional de- 
rivative gives more insights about a disease under consideration 
[11,26,29,30]. Different fractional operators with singular and non- 
singular kernel were suggested in literature [5,27]. The applications of 
these fractional operators can be found in recent literature and refer- 
ences therein [3,8,9,20,33-35]. Recently, very few COVID-19 models 
based on fractional order operators are proposed. The dynamics of 
(2019-nCoV) fractional derivative to explore the situation of the 
pandemic in Wuhan is studied by Khan and Atangana [19]. Some other 
work about the mathematical models and its applications to fractional 
operators can be seen in [28,36,21,6,7,23,18].The authors in [28] 
considered a mathematical model for hepatitis C in fractional derivative 
and examined its results. The CoVID-19 model and its dynamical anal- 
ysis using fuzzy Caputo differential equation is considered in [36]. The 
TB dynamics with relapse has been considered in [21] where the authors 
explored the results using conformable derivative. Some interesting re- 
sults regarding fractional calculus and its connection to attract many 
physical structure is studied in [6]. The authors in [18] obtained the 
stability and numerical results for a HIV/AIDS model in fractional 
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derivative. The oxygen diffusion problem in various fractional operators 
is discussed in [23]. A novel fractional-fractal operator with the singular 
and nonsingular kernel is applied to obtain a dynamical model for 
coronavirus is studied in [4]. Baleanu et al. [10] examined the dynam- 
ical investigations of coronaries model in fractional derivative using 
exponential kernel. A fractional order COVID-19 model in Caputo sense 
is recently studied in [31]. A more recent work on coronavirus consid- 
ering the infected cases observed in Kingdom of Saudi Arabia is inves- 
tigated in [2]. 

Motivated by the above discussion, currently in this work, we study 
the dynamical analysis of the COVID-19 model presented in [32], 
considering the Caputo operator in order to gain more insights about the 
pandemic. We consider here the real infected cases of novel coronavirus 
in Pakistan and obtain the results. A detailed theoretical analysis of the 
fractional order COVID-19 model is presented. Further, we estimate 
model parameters and the basic reproduction number using confirmed 
COVID-19 cases reported from the beginning of the outbreak till June 
30, 2020 in Pakistan. The impact of important model parameters is 
shown for various values of the arbitrary fractional order of the Caputo 
operator. The rest of paper organization is as follows: The basics con- 
cepts related to fractional order derivatives are presented in Section 2. 
The model formulation in integer with parameter estimation and curve 
fitting to reported cases and its generalization to the fractional order 
derivative are presented in Section 3. Section 4 contains the analysis of 
the model and its stability results. In Section 5, we determine the 
important parameters and there impact on the basic reproduction 
number through sensitivity analysis is presented. The graphical results 
through biological discussion is depicted in Section 6. In Section 7, the 
work has been summarized with important suggestions about the dis- 
ease control in the society. 


2. Preliminaries on fractional derivative 


In this section, we briefly discuss background material from frac- 
tional calculus. 


Definition 1. Consider u € C" be a function then the Caputo definition 
with a in (r —1,r) where r € N [27] is defined as: 


Df) =~ f “Ra 


(t — a) 
obviously, “D%(u(t)) approaches to u (t) whenever a1. 


W, 


Definition 2. The fractional integral having of order a > 0 of the 
function u : R* 2 R is 


1 [ um) — m) 5 


iu) = cs], ES 


where 0 «a « 1, t>0. 
Definition 3. A generalized Mittag-Leffler function has been devel- 
oped by Atangana and Baleanu defined as [5]: 


ABC DS (u(t)) = Zo f u (m)E,| — am du, 


where a € [0,1). 


Definition 4. The associated fractional integral of the ABC derivative 
is defined as; 
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Table 1 
Definitions of the variables of the system (1). 
State variable Definition 
S Susceptible people 
E Exposed people 
I Symptomatically-infectious people 
Tq Asymptomatically-infectious people 
Q Quarantined people 
Th Hospitalized people 
I, Critically infected people 
R Recovered people 
N Total population 


Table 2 
Estimated and fitted values along with description of the model parameters. 


Symbol Definition Value/day Ref 

^ Requitement rate p*N(0) Estimated 

r Natural mortality rate 1 [1] 

(365 x 67.7) 

& COVID-19 induced mortality rate in 0.022 [17] 
symptomatic class 

y Transmissibility rate relative to Iq class 0.5620 Fitted 

w Incubation period 0.2691 Fitted 

0 Proportion of the symptomatic infection 0.2350 Fitted 

Y Rate at which symptomatic class are 0.4507 Fitted 
hospitalized 

Ëz COVID-19 induced mortality rate in 0.0131 Fitted 
hospitalized class 

Ëz COVID-19 induced death in critically- 0.1622 Fitted 
infected class 

hı Recovery or removal rate of quarantined 0.4870 Fitted 
class 

hz Recovery or removal rate of hospitalized 0.5431 Fitted 
class 

$3 Recovery or removal rate of critically 0.3838 Fitted 
infected class 

$i Rate of recovery in I 0.6381 Fitted 

95 Rate of Rate of recovery in Iq 0.0803 Fitted 

ó Quarantined rate of exposed individuals 0.3716 Fitted 

B Infection transmission coefficient 0.7806 Fitted 

[4 Hospitalization rate of quarantined 0.1396 Fitted 
individuals 

T Transmissibility of infection relative to Ij 0.4826 Fitted 
class 

h Moving rate of Ij to I, class 0.4593 Fitted 


ABC ya l-a a i aci 
al (u(t)) = ABCA) u(t) + AOC | u(w)(t — e)" dw, 


where a € [0, 1). 


Definition 5. 
y' which is known to be its equilibrium point then, 


[37] Consider a constant point for the Caputo system say 


“Diy(t) = u(ty(r)). a € (0,1), 
if and only if u(t, y") = 0. 
3. The classical integer order model 


In this section, we briefly present description of the COVID-19 model 
in integer case. The model was recently proposed in [32]. To formulate 
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the mathematical model, the whole population at time t, denoted by N(t) 
is further split into eight mutually-exclusive sub-classes including sus- 
ceptible S(t), exposed E(t), those infected people who infected after 
some specific incubation period with symptoms (or symptomatic 
infected) I(t), those have no specific disease symptoms (or asymptoti- 
cally infected individuals) I,(t), quarantined class Q(t), individuals 
which are in hospitals or in self-isolation at home H(t), critically- 
infected or COVID-19 infected people who were in ICU I(t) and 
finally, the individuals successfully recovered from the infection is 
denoted by R(t). In the class I;(t) infected individual showing mild or no 
symptoms and the class Ij(t) contains the patient admitted into the 
hospital and those who are self-isolating at home with clinical symp- 
toms. The state variables of the model are described in Table 1, while the 
complete description of the model parameters are describe in Table 2. 
The flow of among model compartments is depicted in Fig. 1. The 
COVID-19 transmission model rely on these subgroups and at the rate 
they interact, governed the non-linear system of ordinary differential 
equations as follows: 


d$ BU + yla + th) 
d ^ N S — ps, 
dE  f(I ^ yL, + th) 
u N S—(@+6+ NE, 
dl 
g7 OME — (9 +uté +y), 
< t= eon (9 4h, 
(1) 

dQ 
qr 9-06 +90. 

t 
dl, 

g, TYI + EO (td c6 rd), 


dl, 
dt = bly — (u + $5 d &)L; 


dR 
P Vil + Vala + $40 + Poln + ble — UR, 


subjected to initial conditions 


S(0) = So, E(0) 
Q(0) = Qo, 1, (0) 


Eo, 1(0) 
Iro, 1-(0) 


lo, La(0) = lo, 
Lo, R(0) = Ro. 


3.1. Parameter estimations 


The purpose of this section is to estimate the parameters with help of 
the least square curve fitting from the confirmed cases of coronavirus 
that notified in Pakistan since March 1 till 30 June 2020. The recruit- 
ment as well as the natural death rate respectively given by A and y are 
estimated from literature and the rest of the parameter values are fitted 
from the reported data. The detailed estimation procedure can be found 
in [32]. The best fitted model predicted curve is shown graphically in 
Fig. 2 and the model parameters with resulting estimated or fitted nu- 
merical values are shown in the Table 2. The numerical value of the -20 
obtained via the fitted parameters is %o ~ 1.8870, which is in agree- 
ment to the value estimated in recently published article [32]. 


3.2. Model in Caputo sense 
In this subsection, we formulate the fractional order COVID-19 


mathematical model. In order to observe the memory effects, the sys- 
tem (1), in term of integral as: 
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Fig. 1. Flow diagram of COVID-19 model. 


x10 


a Reported cases 
—— Model fit 


Cumulative COVID-19 infected cases in Pakistan 


Time (days) 


Fig. 2. Curve fitting (solid red line) to the confirmed infected cases using model 
(1). The data is consider from 1 March to 30 June 2020. 


dS f. AT. BUS Vid ch) 
E d ^ d 
3r [ (t D| N S — uS dt, 
dE i mI" T, + cl, 1 
Tofu E VL, + th) c (o5 9E df, 
a |, N 
do. 
zf k(t — 1) [OE — (8, +p +E + p EE, 
dt Jv 
d Poo, , 
= k(t — LO — O)@E — (82 + jg), dt , 
dt " 
e (2) 
Be [ kt— Oise - (u+ gi + Nola, 
to 
dl, i ] 1 
q T | Rem tlt + oO - (u+ p +E + PM Jar, 
to 
d. ff, l 
= f e- Dio - (ids edi, 
dt i 
aR off, l 
T= f ke IDI + Bale 0 doli + dle — wR 
to 


where, k(t —t') act as time-dependent kernel and by power-law corre- 
lation function as: 


$ 


t)= 


1 2 


ean, 


(a — 1) 


Substituting the value of kernel in Eq. (2), and then applying the 
Caputo fractional derivative of order a —1, then we have; 


k(t— (3) 


C pei E = patr) |^ Bü ve + Th) g as. 
pe E apie a x Án SH e id 2: ; 
"De H =D UD OeE — (9 +u +é +7), 
“pe! E =DD — 9oE — (85 + Wh), 
dm (4) 
"pe? E = “DYE — (u$ 00], 
"or E] = oirro- eden. 
"pe E] tot eto 
"pr? » = “DET + Sala + PO dol s. — HR, 


Both are the inverse operators then we will get the system (4), 
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BU + wl, + th) 
N 


BI + wl, + th) 
N 


CD's(t) =A S — uS, 


“D°E(t) = 


S—(@+6+p)E, 

CD I(t) = OWE — (9% +u 4- & +71, 

“DT, (t) = (1 — OME — (85 + u)L,, (5) 
"DtQ(r) = 8E — (u + hı + CQ, 

DEI (0) = yl + £Q — (u + b + E + Ol, 

“Dile(t) = bln — (u + ds + Ele, 


CD" R(t) = 91 + 931, + OO + oly + 4I. — HR. 


Let 


— BU + ya + thy) 


a) L—— 


and 


K, = (@ có u), K = (9j +u + +7), Ks =(2 + p), 

K, = (u + hı C) Ks = (u + do + 6 + $), Ks = (u + ds + &). 

Then the novel coronavirus COVID-19 model in Caputo derivative 
form is; 


CD'I(t) = OWE — Kol, 

“DI, (t) = (1 — 0)@E — Kala, (6) 
“Di Q(t) = 8E — K4Q, 

“Dt, (t) = yl + CQ — Ksh, 

“D°I.(t) = bly — Kole, 

©DER(t) = 91+ $41, + P,Q + dI + 4I, — UR 


subjected to the non-negative initial conditions. 
S(0)20, E(0)20, 1(0)20, a(0)20, Q(0)20, 1,(0) 20, 1-(0)20, R(0)20. (7) 


In the above Caputo system (6), the model developed here has eight 
components, S, V, E,L,I;,Q,Ij, I; and R known as COVID-19 model of 
disease transmission also called fractional differential equations (FDEs) 
of eight-dimensional dynamical system. 


3.3. Properties of the model 


We provide in details in the following the model important 
properties. 


3.4. Invariant region and attractivity 


The dynamics of the fractional model (6) will be analyzed in the 
following biologically-feasible region, 


EcR', 
such thats 
E= {(S(t),E(t) I(t) a(t), Q0) In (2) L(t) RO) € R$ SO) + E(t) 4-1(0) H(t) 


A 
+O) +) HC) + ROT}. 
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Lemma 3.1. The region given by ECRÉ, is positively invariant for the 
system (6) with the initial conditions in R$ . 


Proof. We have the following for the model (6), 
CpeN(r = ©D*S(t) + £D'E(t) + £D'I(r) + &D?I,(r) + £D*Q(t) + £D? (t) 
+°D*1,(t) + £D*R(r). 


Hence, 


“DIN (t) = A — uN(r) — I(t) — Eola (t) — Esle(t)<A — uN(), 


clearly we have, 
CD*N(t) -- uN(t)&A. 


By applying the Laplace transform, 


L[CD2N() + uN(0)«LIN], 


A 
s*N(s) — s"! N(0) + uN(s)&-— 
^ gt 
S)< N(O : 
NGST +4) S ^s tu 


Applying Laplace inverse, we obtain, 


N(0)&N(0)E, ( = ut^) oF At Egat ( ie ut"), 


The Mittag-Leffler function describe by infinite power series i.e; 


and laplace transform of Mittag-leffler function is 


gt PB 


L|- ] Eup (£at’)| = 


stra 


So N(t) converges for tco and for all t > 0 the solutions of the model 
with initial conditions in € remain in E. Thus, the region E is positively 


invariant and attracts all solutions in R$. To show that the model has 


positive solution, we consider 


RY = {y € R*|y20] and y(t) = (Sr) E(r), 10), a(t), OC), An(0) L(t), R). 


Corollary 1. [24] Suppose that g(t) € C[a, b] and CDfg(t) € (a, b,where 
a € (0,1]. Then if 


(i)-£D*g(t)20, V y € (a,b), then g(t) is non — decreasing. 


(ii).£D*g(r) € 0, V y € (a, b), then g(t) is non — increasing. 


3.5. Positivity and boundedness 


Proposition 3.1. The solution of the system (6) is non-negative, 
bounded for all (S(0), E(0),I(0), Ia(0), Q(0), I} (0), 1; (0), R(0)) € R$, and 
also defined for all t > 0. 


Proof. In order to explore the solution non-negativity, it is require to 
prove that on each hyperplane bounding the positive orthant, the vector 
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field point R$. Form the system (6) we obtained: 


CDSs = A20, €CD'E|, , = 4920 
CD'l,,-  0eE20, Dill,- = (1 — 0)oE20 
“D Olo- = E20, -D?L|, y = YI  CQ20 
CD, o = i20, 

“DIR|go = HIH 84 + GO + dul + dsL 20. 


Thus, by using the above corollary, the desire goal has been obtained i.e, 
the solution will stay in R$ and so, we can have the following region that 
biologically feasible: 


E — ((S.E,L la; O, Inde, R) € RS :S,E,L I, O, Ig 1c, RZO.) 
Since all terms of the sum are positive, then the solution of system (6) is 
bounded. 


Lemma 3.2. The region E is positively invariant for the model (6) with 
initial conditions in R5. Thus, it is enough to study the dynamics of the model 
in the region given in €. This region can be assumed biologically and epide- 
miologically well-posed for our model considered. 


4. Analysis of the model 
4.1. Disease-free equilibrium DFE 
For the equilibrium points of the fractional model (6), we have 


Cpes(r) “Di Q(t) 


CD E(t) = “D*I(t) 
= “D°R(t) = 0. 


Cper (r) Cper(r) = €DelI.(r) 


The DFE is denoted by E, = (S9, E°, I9, I9, Q°, I, I9, R?) is as follow: 


ES (Ž,0,0,0,0,0,0,0). (8) 


4.2. The basic reproduction number 


To derive the important threshold parameter also known is the basic 
reproduction number ^o, by using the next-generation technique, for 
the dynamics of the disease, whether the infection in population is 


decaying or growing. Let x = (E,I, Ia, Q, In,I¢)", then we have 


dx 
C -g-, 
dt 
where 
"i K\E 
E KI — OME 
P 0 v» | El- (1 — 0)o, 
ida " di K4Q — ôE 
KsI, — yI — Q 
0 Kel. — bln 
0 


The Jacobian of above matrices for the linearized system at disease 
free-state is: 
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0 B fy 0 fc O K, 0 0 0 0 0 
000000 -6o K 0 0 0 0 
fie 000000 ve -(1-0o 0 K 0 0 0 
00 00 0 07’ —ó 0 0 Kk 0 0 
000000 0 -y 0 -č Ks 0 
000000 0 0 0 0 -¢d Ks 


Hence the next-generation matrix is obtained as 


poo | f«(K;oC--Kiy0o) pyoll-0) P Bry py prt pr 


KK; KıK,K;K; KK, K KK; Ky KK; Ks ° 

0 0 0 0 00 

FV"'= 0 0 0 0 00 
0 0 0 0 00 

0 0 0 0 00 

0 0 0 0 00 


The basic reproduction number .Z/, at DFE is obtained by taking the 
spectral radius of the Next-generation matrix p(FV !) for fractional 
model as follow: 

P[K2K366t + K,0(K30(Ks + ty) + KjKsy(l E 0))] 


Ro = 9 
KV KK KA KS (9) 


4.3. Local stability of DFI 


[ca 


Theorem 4.1. Letm > O andn > 0 are the integers such that gcd(m,n) = 
1 for a = © and M = n, then the DFE of the model in fractional derivative is 
locally asymptotically stable if \arg(A)| > 35; for all roots A of the associated 
characteristic equation, 


det(diag A" 4"" 4M jMe Me] — J(Ey )) = 0. (10) 


Proof. For local stability, the linearized system is the jacobian of sys- 
tem (6) at Ey implies that: 


-pu 0 -B -Pwy 0 -f O0 0 

0 —Ki p py 0 pt 0 0 

0 00 -K 0 0 0 0 0 

«| 0 (1-@)o 0  -k; 0 0 0 0 
J(Eo) = 0 ô 0 0 -K, 0 0 0 
0 0 Y 0 ¢ —-KS 0 0 

0 0 0 0 0 $ -K 0 

0 0 8, 8; Qı $» ps -u 


After the evaluation of determinant Eq. (10) implies, 
(am + uy) u” + Kg) n +a" +a" +a," -F a4 A" +as) == 0, (11) 


here we can say that the three eigen values are —j, —4, —Kg have 
negative or negative real part and the co-efficient mentioned above is of 
the form 
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ai= Ki Ky Ky Ki4- K5, 
az= K,K;-- Ku K3- Ko Ks - Ki K4- Ko Ka - KaK4-- Ki Ks -- Ko Ks -- KS Ks 
+KiKs — Ba(03- w(1—0)), 
az= Ky Ko Ks - Ky Ko Ka - Ky KaK4- K KaK4-- Ky KoKs-- Ki K3Ks + K K3Ks 
HK V K4K; - Ka Ka Ks - Ka Ki Ks — foo (y8 -- Ks +K4 + K5) — proc 
—fyo(1— 0)(&» 4- Ki - K5), 
d4— Ky, KoK4Kí4- Ky K2K3Ks5 + K Ko KaKs-- Ki K3K4Ks + Ky K3K4Ks 
—pOwyt(K3-+K4) — Brd¢(Ky+K3) - Byax(1 — (Ko Ks -- KoKs +K4K;), 
üs = Ki K,K3K4Ks — B|OC Ko Ka1 4- Kap(K50(Ks +yt) +KiKsy(1 = 0))], 
= KK KsK,Ks(1 — 2). 


Clearly, a; are all positive if Zo <1. The argument of the root of 
equation (4" + 4)? = 0 and (4" +K6) = 0 are similar, that is: 


m 2m n m 
A) = k—>—> 
args) m * m M 2M 


, where k =0,1,2,....,(m—1). 


In similar fashion, we find out the arguments of the roots of equation 
(49P +a, 49" + ay 49" a5 1?" +a44™ +a5) = O are all greater than zy if 
Zo < 1, having an argument less than 3 for Zo > 1. Thus the DFE is 
locally asymptotically stable for Zo < 1. 


Lemma 4.1. 
Bo >l. 


The disease-free equilibrium of the system (6) is unstable if 


4.4. Global stability of DF! 


[s 


To show global asymptotic stability (GAS) of DFE E of the model (6) 
by Layapunov function approach, we proceed as in the following result. 


Theorem 4.2. 
(8) is GAS. 


If Zo<1 then DFE of the Caputo COVID-19 model given by 


Proof. Let the Lyapunov function be in the form below: 


M RoKs Ks + ty Ksy 4 
P(E, I, la, Q, Ip) = E I T; In. 
F(E, L, la, Q, In) ( pr ) +( Ex ) (= + Ki QI, 


The Caputo fractional derivative of (E, I, Ia, Q, In), along model (6) 
satisfies: 
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RoKs\ c K: > K. 
“DIN Esl dasOss)=( ; Joe ( EET epe. (=) “Dia 


pt 2 3 
ES [4 Cpt Cpe 
K Q+ Di, 
_ Ks | PSU+wl,+tl) m Ks+ty &E—K;I) 
pr N Kot 
Ksy č 
+ E-—K;l, óE—K. 1 
(RE) a - meg (E) (a£ ker 
TCQ—Ksl,, 
RoKs 


< 


Ks+ty 
js B - y c) Kiel ( Kat ) (oe K1) 


(=) (a 0)E Ku) E) (se K,Q)+71 


TCO-K;I,, SSN. 


Ks, K. RK K. 
( "y ‘1)+( 7 vl, yl, Gifs - Kil) 


T T 


Ks+ty 
Kt 


E+( 


66 BoKsK 
(x oKs Kı 


wKs 
OWE. 1—0)o0E 
K, pr jme e jo ) 


K K 
(Fo 114 Po yl, 4-Ks( 225 —1)1;, 


K. 
«(40-1 HL, Fl ), 


<0. 


Hence, it follows that “Diy <0, for Zp <1 all the parameters and 


variables are non-negative with “Des 0 iff E =I =la =Q =I, =0. 


So the infected compartments (E, I, Ia, Q, Iņ) approached zero whenever 
t>o. Using E =I = I; = Q = I, = 0 in equations of system (6), we have 


S^ 1,250 and R0 as too. Thus using lypaunov stability theorems for 
fractional case developed in [37], every solution of the Caputo model (6) 
with non-negative initial conditions, approaches to Ej as too in E. 
Thus, it follows that the DFE of the model (6) is GAS. In the next sub- 
sections, we have to explore the positive equilibria of the model where at 
least one of the infected component is non-zero. 


4.5. Existence of endemic equilibrium point 
The arbitrary endemic equilibrium point (EEP) for the Caputo model 
(6), represented by 


Lon _ Serre o” r R“) 


stea 


with N" —(S" +E" +I" «I, +Q” +I, +I; +R”). Nowat steady 
state, after solving the equations of Caputo model gives 
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ii ^ 
S = ~, 
A+u 
T AA AS" 
E — — 


Kia +p) K?’ 


fe ON 0oE" bwa S* 


I = 7 ais”, 
Ki,KY(A +p) Ky K, ky 
r= a- Ooi A |.(1-60)eE^ (1-0)o/S" - às, 
Ki K3(A +p) K; KiK; 
" 5A A ôE” és" Em 
Q = 7 = = —d4AS, 
KiKi(A +u) K KiKa 
2o gat E“ E 
D= 2 >e L2 (Kym + Kó = (Kyy0w + K-t) 
à mo ag, CEU E EIE 
2 o; ora vs” jas” 
[= == = d(Kay0o + Kd = 
c K, TI XX P(Kay0o + Ka56) e KKK; Ke ` 
T Lor us m m m ay dsA S” 
R” = D(AIT + Bal’ + 10" + daly + dsl) = — 
(12) 
Further, we define the force of infection at endemic state as below: 
» PI +y vc, 
À $ tyl, tuy) (13) 


N“ 
substituting (12) in (13), we have 


_ PORE + y Caer” Qui S- + 2(KyyOw + KC) L3) 


after simplification, we get 


(+N =AZ», 


Hence, 
» a-l ] 
A= >0, >l, 
ae , 0 ; 
where 
1 dd, 
dg = —-F di di + d3 + dy += ++. 
6 [d i doc did +E +5 


Lemma 4.2. There exists a unique endemic state for the model (6), 
whenever 2o > 1. 


4.6. Local stability of (EEP) 
To formulate the Jacobian of system (6) at endemic state. The 


following result is obtained as: 


Theorem 4.3. The endemic equilibrium E," of Caputo-fractional model 
(6) along with N = N"' is (LAS), if Ho > 1. 


Proof. The Jacobian of system (6) at associated endemic equilibrium 


(EE) follows that; 
J(E, ) = (UassVasa]; 


where, 
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se r Suy BS" Sat Bus" sta 


an ON wY wee (Uy) 
stag s) ps” sa" pws sa" 
(Py NT Ki WN ONY N QUY 
0 0o —K, 0 
im 0 (1-0) 0 -K; i 
0 6 0 0 
0 0 
0 0 0 
0 0 vi $5 
and 
S" prs S"4 S"4 S" 
(Ny N^ (NP (N"y | (NY? 
S"A4 pis S" sy sey 
(7 N? (NTY NY (TY 
0 0 -0 0 
= 0 0 0 0 
—K, 0 0 0 
4 -K; 0 0 
0 h — Kg 0 
Qı py 3 -u 


One of the eigenvalue is —, that is negative and others can be obtained 
from the equation: 


A! 4b AS + BAP + b3A* + b44? + bs4? + bgÀ +b; = 0. (14) 


Now the Hurwitz matrices are: 


7 fh 1 
A, =b;, m= (p a 


b 1000 
b 100 
b 10 b, b; b 1 0 
H &d4lhmel lg bs by by by b 
3 E h à , 4 bs by bi b; , 3 om Du ‘ a | 
5 V4 03 0 0 bs b, 5 U4 03 
000 0 bs 
b 100000 
b 10000 
b b b 1 0 0 0 
b, b b 1 0 0 
bs by by by bi 0 
bs by by by by 1 
He = .H;—| 0 0 b, by by by bi 
2-0 Ao HS de 0 0 0 0 bs by b 
0 0 0 0 bs b, Mb 
0 0 0 0 0 bs bs 
0 0 0 0 0 b 
00000 0 b 


Here we can say that the coefficient signs b; are strictly positive for the 
roots of the polynomial to be negative. All roots of polynomial have 
negative real part iff, DetH; > 0, j = 1, 2, ...., 7. Thus, the conditions 
ensures the endemic stability. 
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The basic reproduction number, 


(a) (b) 


Fig. 3. Impact of effective contact rate / and quarantine rate ô on “po with contour plot. 


4.7. Global stability of EEP 
A^ = BST, 
For global stability at endemic state of the Caputo COVID-19 model KE" = AS", 

(6) is given for the special case when infected hospitalized do not bl" = aug" (18) 
transmit infection (t = 0), and disease-induced mortality is negligible Eg = do 6o EU 
(& = é = 0), then the reduced model is obtained as: Ya 7 

Cpes() =A pus We uc Y X =p ty, &e, duy 

Cpag(o) — PU E Wie Fh) S to 4 54 pE, 


N 
CD'I(t) = OWE — (81 +u +y), 


Theorem 4.4. The fractional model (15) at unique EEP ie., E, * * is GAS 
in E, if 91 > 1 and the condition (25) holds. 


Proof. Consider the following non-linear Lyapunov function for the 


Cpa = EB = 

FG DOE NS (15) model (15), so that the corresponding unique EEP E; ' * exist by letting 
“Di Q(t) = 8E — (u + d, + OQ, Ry 7 1. 

“Dth(t) = y - CQ — (u+ drt bh, R(t)= (so S sn) (eo E" Ene) 

Cha = - 

DIT -(t) = bln — (p + ds + £), Ki (ro rere O). K, (: ()-r"-£"1 ae) 
CDIR(t) = 911 + Bala + 0 + bol + dale — HR, o 7) (Lo t t nU 


where, now 
The time fractional derivative of (19) is: 


A(t) = BU + Wha) . 
^ CDAR(t) = (1 B pesto (1 B Jeore es (1-1) Epere) 
The reduced model shows that N(t)>fas too. Finally we have 5 i o L 
Ki D ex 
Al) = Ath), p=". (16) ‘apa! I, ) DA) 
The relative expression for the reproduction number o, of the 7 SU E" 
reduced model with (16) is (: S J^ 5-4 (: E Jes KE) 
„ _ Bio|K30 + ky(1 — 6)] Kf, r7 K p 
Zu = € 17 1 1 a 
01 [qw (17) at j ) (awe kol)4 me L Ja 0)@E—K3l,). 


The following results for model (15) at steady state obtained as: 
Using the solution from (18) gives 
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Fig. 4. Impact of effective contact rate // and hospitalization/self-isolation rate y on #%o with contour plot. 


The basic reproduction number, 
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Fig. 5. Impact of transmission coefficient y and quarantined rate ô on the basic reproduction number 7, with contour plot. 


SN S hee ee T" 
1 D'S(r)— (1 I--I,)S—pS), E E pb ghee dx avus TE 
jJ ug Bro), — (, E eee (1 PES-A S ev $7) E 
E E E 
=p,S T" sus E  E"'IS IS 
Bis n =pS I 1—— LE ba 
E" Egg" [Us 
— S" LS Ll m 
Pys DL, Ll——L—ue x UE: E E"LS LS 
4 S I I HYSTE” | l—--————ás mm. 
«5 da PWS tq E" ELS" rs 
"T S" sS (21) 
4S (2————— ]. 
à S s 
(20) 


10 


A. Ali et al. Results in Physics 20 (2021) 103669 


The basic reproduction number,R, 
1 
03 
b 05 


14 

= 

, 10 
0 01 02 03 04 05 06 07 08 08 1 


04 
03 
02 
01 

0 


(a) (b) 


Fig. 6. Effect of transmission coefficient 7 and quarantined rate ô on the basic reproduction number .#o with corresponding contour plot. 
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Fig. 7. Effect of hospitalization/self-isolation y and quarantined rate ô on the basic reproduction number .Z with corresponding contour plot. 
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Time(days) 


(a) 


Time(days) 


(c) 


Substitute (20) to (23) in (20), after simplification gives 


"ors As n (a ELLE EE) 
S L" ES" IE" LE" E 
E"LS IE LE 


=) 
El, S IE LE E 


(24) 


Since the arithmetic mean exceeds the geometric mean, we have the 
following interpretation form (24): 


j--6 sy 


Stt S 
2 TEES 
( S St 


Further, if the following inequalities holds 


<0. 


SUL jue. gem gm E 
4 oe ee ae — oe tir | <0 
S Ln ES IE LE 
x E - N (25) 
S^ I E"LS I"E I"E E 
4 = mE ae — oe tion | <0 
S I EL” S IE LE E 


then °D?R(t) <0 for #1. Thus, &(t) is a Lyaponov function in E. 
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Fig. 8. Dynamics of fractional COVID-19 model for different values of a. 


Therefore, by Lasalle invariance principle, that 
ims = 
=i,’ timo = QU 
2n 


t limE( =E "", lim/(t) = **, limL(f) 
t>o too too 


^, limh) =F, ^, limL(t) 
100 lim 
E limR(t) —R** 
t>00 


Thus every solution of the reduced model, tends to its unique endemic 
equilibrium for associated reproduction number as t— co. For the GAS of 


EEP (Ej ) whenever #o > 1, through the conjecture. 


Remark 1. The unique EEP of system (6) is GAS in E, if Hp > 1, in E. 


5. Zo versus model parameters 


Figs. 3-7 depict the impact of model various important parameters 
on the value of .#o with the corresponding contour plots. The variation 
in the value of Zo versus the effective contact rate / and quarantine rate 
6 is shown in Fig. 3. It is observed that .Z/ decreases for values smaller 
than 1 with a decrease in J and an increase in ô. Similarly, Ho can be 
reduced to a smaller value less than 1 with the in 
hospitalization/self-isolation rate y and decrease in the effective contact 
rate f. This interpretation can be found in Fig. 4. The combined effect of 
the infectious rate due to hospitalized infected patients y and 6 is 
analyzed in Fig. 5 while the behavior of .Zo for the variation in t and ô is 
depicted in Fig. 6. This interpretation demonstrates that .Z can be 


increase 
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Time(days) 


(a) 


Time(days) 


(c) 


effectively reduced with a reduction in the effective infectious rates (7 
and ó) and enhancing the efficacy of contact-tracing and quarantine 
policies of the exposed population. Finally, we analyzed the role of 
variation of both quarantine rate y and hospitalization rate y on the 
value of #9. It is worth mentioning that o remains very smaller with 
an increase in both these parameters. 


6. Numerical results and discussion 


This section presents the numerical results of the Caputo COVID-19 
model (6). The generalized predictor-corrector of Adams-Bashforth 
Moulton method developed in [13,11] is considered for the numerical 
solution of the model (6). To simulate the model, we use the real pa- 
rameters fitted given in Table (2). To study the model qualitative 
behavior and its parameters effect on the system and possible control, 
we varied the model important parameters and the fractional order a 
and obtained the simulation results. We studied graphically the coro- 
navirus model (6) for five different values of a are shown in Figs. 8 and 9. 
It is observed that the susceptible population is increased while the 
population in all remaining classes is decreased significantly, when the 
fractional order a takes smaller values. The impact of effective contact 
rate 2 on the dynamics of symptomatic and asymptomatic infected in- 
dividuals is shown in Fig. 10 for four different values of fractional order 
a of the Caputo derivative. A significant decrease in newly reported 
symptomatic and asymptomatic infected individuals is seen with a 
decrease in J. Further, with a 40% decrease in effective contacts (i.e., 
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Fig. 9. Dynamics of fractional COVID-19 model for different values of a. 


pf = 0.4683) to its estimated baseline value, a dramatic decrease is 
observed in the peak of infection curves. This trend is observed for all 
values of a, but for smaller values of fractional order, the decrease in the 
infection curves is comparatively more significant. Fig. 11 illustrates the 
influence of effective contacts # on the behavior of hospitalized and 
critically infected (or ICU patients). It is seen that the total hospitalized 
and the critically infected population is decreased very well when f is 
decreased at 1096, 2096, 3096 and 4096 respectively to its baseline value 
as shown in Fig. 11. Moreover, the decrease in the cumulative hospi- 
talized and critically infected individuals becomes more significant for 
smaller values of fractional order a. As a consequence of this graphical 
interpretation, the severity of COVID-19 infection can be reduced by 
following a strict Standard Operating Procedure (SOPs) in order to 
reduce the effective contacts among infected and susceptible in- 
dividuals. The impact of quarantine rate or contact-tracing parameter ó 
on the dynamics of symptomatic and asymptomatic infected individuals 
is analyzed in Fig. 12. From this analysis one can observe that with an 
increase in 6 (up to 50 96 i.e. 6 = 0.5574) to the baseline value a sharp 
decrease in the peak of infected curves is seen as shown in Fig. 12. We 
depicted the simulations for four distinct values of a and obtained 
identical behavior for fractional order parameter values, that the 
decrease in the infected individuals is comparatively faster for smaller 
values of a. Furthermore, Fig. 13 depicts the dynamics of the cumulative 
hospitalized and critically infected population for variation in the 6. 
With the increase in the contact-tracing policy, a decrease in hospital- 
ized and the critically infected population is observed. The same 
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Fig. 10. Influence of ó (quarantine rate) on the symptomatic and asymptomatic infected individuals with different value of a. 
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Fig. 11. Impact of f (contact rate) on the hospitalized and ICU infected individuals with different value of a. 
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12. Influence of ô (quarantine rate) on the symptomatic and asymptomatic infected individuals with different value of a. 
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g. 13. Influence of ó (quarantine rate) on the hospitalized and ICU infected individuals with different value of a. 
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patterns are found for all values of fractional order a. Finally, we 
analyzed the influence of hospitalization of confirmed COVID-19 cases 
(y on the new cumulative symptomatic individuals. The resulting 
graphical interpretation for four distinct values of a is shown in Fig. 14. 
A small variation in the proportion of symptomatic individuals occurs 
even if we increase y at a significant rate. It is found that when a de- 
creases, the variation in the infected population is slightly feasible. 


7. Conclusion 


Fractional order derivatives are more prominent and provide 
comparatively deeper insights about real-world problems. Moreover, 
the derivative with fractional order is capable to capture the fading 
memory and crossover behavior found in most biological phenomena 
including infectious diseases. In this paper, we studied the dynamics of 
the COVID-19 outbreak in Pakistan using a fractional order mathemat- 
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Fig. 14. Influence of y (isolation/hospitalization rate) on symptomatic COVID-19 infected individuals with different value of a. 


ical model in Caputo sense which is the most commonly used operator in 
the modeling approach. Initially, the Caputo fractional COVID-19 model 
is formulated and then provided its basic and necessary mathematical 
properties that include the existence and positivity of the system and its 
solution. The detailed stability results both local and global for DFE and 
EE are explored using fractional stability approaches. The impact of 
variation in various model important parameters on the dynamics of the 
basic reproduction number is depicted graphically. Furthermore, the 
model parameters and the basic reproduction number -#9 are estimated 
from the COVID-19 real infective reported cases from health ministry of 
Pakistan from the start of the outbreak till June 30, 2020. The non-linear 
least square procedure is used for the estimation purpose. The present 
finding showed that the results of predictions obtained through our 
model is inline with those published in literature and thus more 
appropriate values of the parameters are estimated. The numerical value 
of the basic reproduction number evaluated using the parameters 
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obtained from fitting is Zo ~ 1.8870 which is inline to the published 
result in [32]. Finally, we simulated the Caputo model in order to 
examine the impact of various NPIs on the dynamics and control of 
COVID-19 in Pakistan. We varied the relevant parameters at different 
levels to its baseline values and obtained the graphical results for 
distinct values of fractional order a of the Caputo derivative. It is found 
that the reduction in the effective contacts / (up to 40 96) and an 
enhancement in contact-tracing policy to quarantine the exposed in- 
dividuals ó (up to 5096) to its baseline value that dramatically reduced 
the peak of infected curves. Moreover, the reduction in the infected 
population becomes more significant for smaller values of fractional 
order a. Thus we conclude that the results obtained for the fractional 
case are reliable, realistic and more biologically feasible. In future, the 
present study can be extended to fractional order with the non-singular 
kernel. 
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